Prospecção Séries Temporais

G4: Edinaldo de Alencar / Igor Freire / Ramon Araújo / Ricardo Ribeiro

30 de outubro de 2014

Importação da base de dados

Nota: a informação sobre os dias da semana é redundante, pois pode ser obtida através de:

dia <- weekdays(index(dataset_diario))

Por isso, na leitura da base em dataset_diario, esta coluna foi ignorada.

Séries Temporais

Conversão dos dados para o formato de série temporal no R.

Série Mensal:

# Frequency --> número de observações por unidade de tempo 
# define a unidade de tempo (e.g. 12: unidade de tempo = ano)
tsMensal <- ts(dataset_mensal, frequency=12, start=1992)

plot of chunk unnamed-chunk-6

Série Diária:

Série temporal diária

tsDiario <- msts(dataset_diario, start=c(2002,1,1), seasonal.periods=c(7, 365.25))

plot of chunk unnamed-chunk-8

Nota: é importante observar os anos bissextos, 2004 e 2008.

Inspeção dos Dados

Série para um dia específico da semana

plot of chunk unnamed-chunk-9

Série para um dia e mês específico da semana

plot of chunk unnamed-chunk-10

Observar a influência dos feriados. Exemplo, 12 de outubro foi uma quinta-feira em 2006.

Inspeção dos Dados - 2

A série temporal pode ser decomposta em:

Fluxo Mensal

plot of chunk unnamed-chunk-11

Fluxo Diário

plot of chunk unnamed-chunk-12

Observações

  1. Clara tendência (trend) de subida.
  2. Sazonalidade de 1 ano.
  3. No fluxo diário, há sazonalidade semanal.

Inspeção dos Dados - 3

Fluxo Diário para um intervalo arbitrário de 8 semanas

plot of chunk unnamed-chunk-13

Inspeção dos Dados - 4

Auto-correlação

acf(dataset_diario, col="red")

plot of chunk unnamed-chunk-14

Parcial:

pacf(dataset_diario, col="red")

plot of chunk unnamed-chunk-15

Objetivos

Realizar prospecções de curto, médio e longo prazo.

Curto Prazo

Para curto prazo, será utilizada a série temporal com dados diários (tsDiario).

Médio Prazo

Longo Prazo

Para médio e longo prazo, será utilizada a série temporal com dados mensais (tsMensal).

Separação dos dados

Separação das séries temporais em conjuntos de treinamento e de testes.

Curto Prazo (tsDiario)

Conjunto de Treinamento

tsDiarioTrain <- window(tsDiario, end=c(2008,183)) # até 30/06/2008

plot of chunk unnamed-chunk-17

Conjuntos de Teste

30 dias:
# De 01/07/2008 a 30/07/2008
tsDiarioTest30Days <- window(tsDiario, start=c(2008,184), end=c(2008,213)) 
45 dias:
# 45 dias a partir De 1/1/2008
tsDiarioTest45Days <- window(tsDiario, start=c(2008,184), end=c(2008,228)) 

Médio/Longo Prazo (tsMensal)

Conjunto de Treinamento

tsMensalTrain <- window(tsMensal, end=c(2007, 12)) # De Jan 1992 a Dez 2007

plot of chunk unnamed-chunk-24

Conjuntos de Teste

4 meses:
tsMensalTest4mth <- window(tsMensal, start=2008, end=c(2008, 4)) # De Jan 2006 a Mar 2006
6 meses:
tsMensalTest6mth <- window(tsMensal, start=2008, end=c(2008, 6)) # De Jan 2006 a Jun 2006
1 ano:
tsMensalTest1yr <- window(tsMensal, start=2008, end=c(2008, 12)) # De Jan 2006 a Jan 2007
2 anos:
tsMensalTest2yr <- window(tsMensal, start=2008, end=c(2009, 12)) # De Jan 2006 a Jan 2008

Bibliotecas Utilizada

Prospecções Estatísticas

ver (R. J. Hyndman 2014) e (Leek 2014):

library(forecast)

Redes Neurais:

library(nnet)

Estratégia de Avaliação dos Métodos

Apresenta-se a seguir:

Nota:

Acurácia

MAPE:

Diferentemente do ME, MAE e RMSE, o medidor MAPE é independente da escala.

\[M = \frac{1}{n} \sum \limits_{t=1}^{n} \left| \frac{A_t - F_t}{A_t} \right|,\]

onde \(A_t\) é o valor verdadeiro, \(F_t\) é o valor predito e \(n\) é o número de amostras temporais preditas.

Como desvantagem, o critério MAPE tem comportamento indefinido quando o valor verdadeiro \(A_t=0\) (divisão por zero).

Benchmark Métodos Estatísticos

Curto prazo

30 dias

##           Training set Test set
## Mean            14.656   22.506
## Näive            6.214    4.674
## N-Sazonal        8.720    9.066
## Drift            6.216    4.788
## Regressao        5.868    7.383

Médio prazo (4 meses)

##           Training set Test set
## Mean             25.58   35.295
## Näive             4.44    5.639
## N-Sazonal         7.68    6.051
## Drift             4.47    6.542

Longo prazo

1 ano

##           Training set Test set
## Mean            25.582   39.867
## Näive            4.440    5.708
## N-Sazonal        7.680    6.898
## Drift            4.470    4.722
## Regressao        4.741   10.023
## ETS              1.442    2.603

2 anos

##           Training set Test set
## Mean            25.582   41.232
## Näive            4.440    6.801
## N-Sazonal        7.680    9.046
## Drift            4.470    5.210
## Regressao        4.741   10.398
## ETS              1.442    2.964

Modelo auto-regressivo (AR)

Em um modelo de auto-regressão, prevemos a variável de interesse usando uma combinação linear dos valores passados da variável.

Assim, um modelo auto-regressivo de ordem \(p\), denotado por \(AR(p)\), pode ser escrito como:

\(y_t = c + \varphi_1 y_{t-1} + \varphi_2 y_{t-2} + \dots + \varphi_p y_{t-p} + e_t,\)

onde \(c\) é uma constante, \(e_t\) o ruído branco, e, os valores defasados de \(y_t\) são os preditores.

Para um modelo \(AR(1)\), tem-se \(\varphi_1\in(-1,1)\)

Para um modelo \(AR(2)\), tem-se \(\varphi_1,\varphi_2\in(-1,1)\) com \(\varphi_1 + \varphi_2 < 1\) e \(\varphi_2 - \varphi_1 < 1\)

Modelos médias-móveis (MA)

Ao invés de usar valores passados da variável de previsão em uma regressão, o modelo de médias-móveis usa erros de previsão do passado em um modelo de regressão-like.

\(y_t = c + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} + \dots + \theta_q e_{t-q}.\)

Esse modelo tem ordem \(q\) e é denotado por \(MA(q)\).

Para um modelo \(MA(1)\), tem-se \(\theta_1\in(-1,1)\)

Para um modelo \(MA(2)\), tem-se \(\theta_1,\theta_2\in(-1,1)\) com \(\theta_1 + \theta_2 > -1\) e \(\theta_1 - \theta_2 < 1\)

Relação entre modelos AR e MA

É possível escrever o modelo \(AR(p)\) como um modelo \(MA(\infty)\). Em particular para um modelo \(AR(1)\), vem:

\(y_t = \varphi_1 y_{t-1} + e_t\)

\(y_t = \varphi_1 ( \varphi_1 y_{t-2} + e_{t-1} ) + e_t\)

\(y_t = \varphi_1^2 y_{t-2} + \varphi_1 e_{t-1} + e_t\)

\(y_t = \varphi_1^3 y_{t-3} + \varphi_1^2 e_{t-2} + \varphi_1 e_{t-1} + e_t\)

e assim por diante. Onde \(\varphi_1\in(-1,1)\)

Modelo ARIMA

ARIMA (Autoregressive Integrated Moving Average)

Equação de Predição

\[y(t) = \sum_{i=1}^{p}\varphi_ix(t-i) + \sum_{j=1}^{q}\theta_j\epsilon(t-j)\]

Notação

\(ARIMA(p,d,q)\), onde:

Método ARIMA em prospecções de médio prazo

Para 4 meses

mesesAPrever <- 4

Modelo:

arima_model_mensal <- auto.arima(tsMensalTrain)

Predição:

fcast_arima_4mth <- forecast(arima_model_mensal, h=mesesAPrever)
plot(fcast_arima_4mth, xlab="Anos", ylab="Fluxo Mensal", include=36)
lines(tsMensalTest4mth, col="red")

plot of chunk unnamed-chunk-41

Acurácia:

accuracy(fcast_arima_4mth, tsMensalTest4mth)
##                    ME  RMSE   MAE    MPE  MAPE   MASE      ACF1 Theil's U
## Training set    716.4 11391  8353  0.115 1.625 0.2089 -0.001948        NA
## Test set     -12192.3 18058 15355 -1.513 1.933 0.3840 -0.632429    0.2873

Para 6 meses

mesesAPrever <- 6

Predição:

fcast_arima_6mth <- forecast(arima_model_mensal, h=mesesAPrever)
plot(fcast_arima_6mth, xlab="Anos", ylab="Fluxo Mensal", include=36)
lines(tsMensalTest6mth, col="red")

plot of chunk unnamed-chunk-45

Acurácia:

accuracy(fcast_arima_6mth, tsMensalTest6mth)
##                    ME  RMSE   MAE    MPE  MAPE   MASE      ACF1 Theil's U
## Training set    716.4 11391  8353  0.115 1.625 0.2089 -0.001948        NA
## Test set     -10206.2 15249 12315 -1.258 1.539 0.3080 -0.585487    0.2658

Método ARIMA em prospecções de longo prazo

Para 1 ano

mesesAPrever <- 12

Predição:

fcast_arima_1yr <- forecast(arima_model_mensal, h=mesesAPrever)
plot(fcast_arima_1yr, include=48)
lines(tsMensalTest1yr, col="red")

plot of chunk unnamed-chunk-49

Acurácia:

accuracy(fcast_arima_1yr, tsMensalTest1yr)
##                  ME  RMSE   MAE    MPE  MAPE   MASE      ACF1 Theil's U
## Training set  716.4 11391  8353 0.1150 1.625 0.2089 -0.001948        NA
## Test set     4503.2 18305 15764 0.4428 1.841 0.3943  0.403569    0.4769

Modelo TBATS

Justificativa

Modelo TBATS, a partir da série multi-sazonal do conjunto de treinamento:

tbats_model <- tbats(tsDiarioTrain)
#tbats_model <- tbats(window(tsDiarioTrain, start=2005))

Decomposição do fluxo diário em termos trend+season pelo método TBATS

plot(tbats_model)

plot of chunk unnamed-chunk-52

O método TBATS é indicado para a previsão do fluxo diário, pois:

  1. Sazonalidades múltiplas: semanal e anual

  2. Sazonalidade semanal é de alta frequência

Método TBATS em prospecções de Curto Prazo

Para 30 dias

diasAPrever <- 30

Predição:

tbats_fc30Days <- forecast(tbats_model, h=diasAPrever)
plot(tbats_fc30Days, include = 120)
lines(tsDiarioTest30Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-54

Acurácia:

accuracy(tbats_fc30Days, tsDiarioTest30Days)
##                  ME  RMSE   MAE     MPE  MAPE   MASE       ACF1 Theil's U
## Training set -13.13 742.6 436.1 -0.1878 2.095 0.2186 -0.0007912        NA
## Test set     412.21 689.8 601.6  1.3682 2.122 0.3015  0.6355962    0.3331

Método TBATS em prospecções de Curto Prazo

Para 45 dias

diasAPrever <- 45

Predição:

tbats_fc45Days <- forecast(tbats_model, h=diasAPrever)
plot(tbats_fc45Days, include = 120)
lines(tsDiarioTest45Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-57

Acurácia:

accuracy(tbats_fc45Days, tsDiarioTest45Days)
##                  ME  RMSE   MAE     MPE  MAPE   MASE       ACF1 Theil's U
## Training set -13.13 742.6 436.1 -0.1878 2.095 0.2186 -0.0007912        NA
## Test set     572.57 869.5 732.4  1.8836 2.518 0.3670  0.6215786    0.3952

Exponential Smoothing

Exponential Smoothing - 2

Ver (R. J. Hyndman and Athanasopoulos 2013). O autor possui um livro exclusivamente sobre o assunto em (R. Hyndman et al. 2008).

ETS

Para cada um dos 15 métodos, há dois modelos possíveis: com erro aditivo e com erro multiplicativo.

Há, portanto, 30 métodos distintos.

Exponential Smoothing no R

Método Exponential Smoothing em Médio Prazo

Para 6 meses

mesesAPrever <- 6

Modelo:

etsMensal <- ets(tsMensalTrain)

Prospecção:

fcastMensal6mth <- forecast(etsMensal, h=mesesAPrever)

Comparar predição com valores reais do conjunto de teste:

plot(fcastMensal6mth, xlab="Anos", ylab="Fluxo Mensal", include = 36);
lines(tsMensalTest6mth, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-62

Acurácia:

accuracy(fcastMensal6mth, tsMensalTest6mth)
##                  ME  RMSE  MAE      MPE  MAPE   MASE     ACF1 Theil's U
## Training set  856.7  9849 7308  0.13394 1.442 0.1828 -0.02704        NA
## Test set     -755.7 10937 8653 -0.06425 1.104 0.2164 -0.63165    0.2848

Prospecção a Longo Prazo

Transformação Box-Cox

Há um fator multiplicativo no fator sazonal. Uma transformação Box-Cox demonstrou ser apropriada.

Série original

plot(tsMensalTrain, xlab="Anos", ylab="Fluxo Mensal")

plot of chunk unnamed-chunk-64

Série transformada

Calcula-se o fator \(\lambda\) para a transformação:

lam <- BoxCox.lambda(tsMensalTrain)

Série Resultante:

tsMensalBoxCox <- BoxCox(x = tsMensalTrain, lam)
plot(tsMensalBoxCox, col="red")

plot of chunk unnamed-chunk-67

Prospecção com decomposição sazonal de Loess - STLF

  1. Aplica-se a decomposição STL
  2. Remove-se a sazonalidade da série temporal
  3. Um modelo é treinado na série resultante
  4. Faz-se a prospecção
  5. Adiciona-se novamente o último período da sazonalidade estimada nos resultados
mesesAPrever <- 24

Dados com sazonalidade removida:

tsMensal.stl <- stl(tsMensalTrain[,1], s.window=12)
# Seasonally adjusted data constructed by removing the seasonal component.
plot(seasadj(tsMensal.stl))

plot of chunk unnamed-chunk-69

stlf_model <- stlf(tsMensalTrain[,1], lambda = lam)
stlf_fcast <- forecast(stlf_model, h=mesesAPrever)
plot(stlf_fcast, xlab="Anos", ylab="Fluxo Mensal")
lines(tsMensalTest2yr, col="red")

plot of chunk unnamed-chunk-70

Acurácia:

accuracy(stlf_fcast, tsMensalTest1yr)
##                   ME  RMSE   MAE      MPE  MAPE   MASE     ACF1 Theil's U
## Training set  -25.42  8186  6175 -0.01584 1.238 0.1544 -0.02348        NA
## Test set     8132.24 17138 13696  0.92760 1.619 0.3425 -0.11417    0.4831

Referências

De Livera, Alysha M, Rob J Hyndman, and Ralph D Snyder. 2011. “Forecasting Time Series with Complex Seasonal Patterns Using Exponential Smoothing.” Journal of the American Statistical Association 106 (496): 1513–1527.

Hyndman, Rob J. 2014. “Forecasting Time Series Using R.” http://robjhyndman.com/talks/MelbourneRUG.pdf.

Hyndman, Rob J, and George Athanasopoulos. 2013. Forecasting: principles and Practice. OTexts.

Hyndman, Rob, Anne B. Koehler, J. Keith Ord, and Ralph D. Snyder. 2008. Forecasting with Exponential Smoothing: The State Space Approach (Springer Series in Statistics). 2008th ed. Springer. http://amazon.com/o/ASIN/3540719164/.

Leek, Jeffrey. 2014. “Practical Machine Learning - Lecture 27 - Forecasting.” Johns Hopkins Bloomberg School of Public Health. https://d396qusza40orc.cloudfront.net/predmachlearn/027forecasting.pdf.